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ABSTRACT 

We address the problem of discretizing continuous cosmological signals such as a galaxy 
distribution for further processing with Fast Fourier techniques. Discretizing, in particular 
representing continuous signals by discrete sets of sample points, introduces an enormous 
loss of information, which has to be understood in detail if one wants to make inference from 
the discretely sampled signal towards actual natural physical quantities. We therefore review 
the mathematics of discretizing signals and the application of Fast Fourier Transforms to 
demonstrate how the interpretation of the processed data can be affected by these procedures. 
It is also a well known fact that any practical sampling method introduces sampling artifacts 
and false information in the form of aliasing. These sampling artifacts, especially aliasing, 
make further processing of the sampled signal difficult. For this reason we introduce a fast and 
efficient supersampling method, frequently applied in 3D computer graphics, to cosmological 
applications such as matter power spectrum estimation. This method consists of two filtering 
steps which allow for a much better approximation of the ideal sampling procedure, while 
at the same time being computationally very efficient. Thus, it provides discretely sampled 
signals which are greately cleaned from aliasing contributions. 

Key words: large scale structure of the universe - data analysis - numerical - statistical 



1 INTRODUCTION 

Many cosmological applications, like estimating the power spec- 
trum from galaxy surveys or calculating the gravitational potential 
of the dark matter distribution in numerical simulations, rely on the 
use of Fast Fourier Transform (FFT) techniques. They are favored 
above other computational techniques as their low computational 
complexity allows to process huge datasets in reasonable times. 

However, since by definition the FFT operates only on finite 
sets of discrete sample points the problem arises how to discretize 
and digitize a continuous signal, in order to record it in computer 
memory and make it available for further processing with FFTs. 
Natural physical signals, though, are generally neither discrete in 
real-space nor in Fourier-space. 

Reducing such a signal to a set of finite and discrete sample 
points will naturally introduce an enormous loss of information. 
For this reason we try to answer what kind and how much informa- 
tion of the true physical system is still represented by the sampled 
data. This is equivalent to ask for a response function describing 
the operation of the computer data acquisition. 

According to Shannon's theorem sampling a continuous sig- 
nal can be achieved by low-pass filtering and sam pling the contin- 
uous signal at discrete positions Jshanno nl ll948lll94g) . Low-pass 
filtering requires to convolve in real-space with the sinus cardinal 
function, which is infinite in extend, and is therefore not practical 
for real world applications. For this reason one usually relies on 



filter approximations of the ideal low-pass filter, which only have 
finite support in real-space and hence allow for fast computation of 
the low-pass filtering convolution. 

In cosmology some frequently used techniques to discretize 
a continuous distribution of point particles are Nearest Grid Point 
(NGP), Cloud In Cell (CIC) or Triangular Shaped Clouds (TSC) 
(Hocknev & Eastwoodll 19881) . which all attribute a weighted frac- 
tion of the individual point particle mass to the surrounding dis- 
crete grid positions. The process of sampling by relying on these 
filter approximation in general not only irrevocably reduces the in- 
formation content of the continuous signal to that one represented 
by a set of discrete points , it also introduces sampling artifacts like 
aliasing or Gibb s ringing dWolberdl997hlMarschner & Lobbll994 
ICui et alj|2008f) . It is often the leakage of aliasing power from the 
stop-band which makes further processing of the sampled signal 
difficult. For example, estimation of higher-order spectra, like the 
bi- or tri-spectrum, with FFT tech niques will be er roneous due to 
sampling artifacts. As mentioned in lCui et ail i2008h . so far, there is 
no known approach to accurately correct for these sampling effects 
in measuring higher-order spectra like the bi-spectrum with FFT 
techniques. The aim for signal processing technologies therefore is 
to find low-pass filter approximations which sufficiently suppress 
these sampling artifacts, while at the same time still being compu- 
tationally far less expensive than the ideal sampling procedure. The 
Digital Signal Processing (DSP) literature provides plenty of pos- 
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sible approaches (see e.g. Smith 2002) which could be introduced 
to various cosmological signal processing problems. 

For cosmological applications recently different approaches 
have been presented to deal wit h the problem of sampling arti- 
facts Jjindl2005l ; ICui et al.ll2008T) . Ijind 12003) proposes to evalu- 
ate all aliasing sums of the NGP, CIC and TSC method. He finds 
that the sampled shot noise contribution to the power spectrum 
can be expressed in a simple analytic way. To correct for further 
aliasing contributions, he introduces an iterative correction method 
for the power spectrum, assuming a power-law behavior. However, 
Jing's method is only applicable to power spectrum estimations 
and in cases where there are no other sources of mode-coupling, 
like masks o f a galaxy surve y or selection functions. The method 
proposed by ICui et al] J2008h utilizes the scale functions of the 
Daubechechies wavelet transformation as filter approximations. It 
has been demonstrated to yield better results than standard CIC or 
TSC methods. 

However, most of these approaches rely on increasing the spa- 
tial support of the filter approximation to more closely represent the 
properties of the ideal low-pass filter. This approach is well known 
in the DSP literature, which favors the windowed sin e functions as 
optimal approximations to the ideal low pass filte r dTheufil et al.l 
l200d : iDuan et alj 120031 ; iMarschner & Lobbl 1 1994 ISmithl 120021) . 
These approaches, though, become computationally more expen- 
sive as the spatial support of the filter approximation grows, and 
tends towards evaluating the ideal low-pass filtering convolution. 
Many applications, like for instance particle mesh codes, however, 
rely on fast but nevertheless accurate low-pass filtering procedures. 
Hence, sampling a continuous signal like the galaxy distribution on 
a discrete grid to further process it via FFT techniques is always a 
trade off between computational efficiency and quality. 

For this reason we introduce a new supersampling technique 
frequently applied as an anti-aliasing technique in 3D computer 
graphi cs to cosmological applications like power spectrum estima- 
tion dWolberdfl^lGoss & Wv]|2000h . 

Our method is a two step filtering process in which the signal 
first is sampled to a discrete grid with super resolution via the CIC 
or TSC method, then low-pass filtered with the ideal discrete low- 
pass filter, and finally resampled at target resolution. 

In this fashion, by using simple low-pass filter approximation 
and utilizing FFTs, we provide a fast method to accurately calculate 
a low-pass filtering procedure. 

As our method uses existing finite support filter approxima- 
tion, it can be understood as complementary to the approach of 
finding better low-pass filter approximations. The combination of 
both approaches will naturally lead to more optimal results. How- 
ever, the supersampling method in itself is computationally much 
less expensive than increasing the spatial support of the filter ap- 
proximation, while at the same time presenting better results as can 
be demonstrated in the case of cosmological matter power spectrum 
estimation. 

Beside being a numerically efficient method, the supersam- 
pling procedure does not only aim at correcting aliasing effects in 
the power spectrum, it also greately cleans the discretized signal 
from aliasing contributions. 

This is important for applications which operate on the sam- 
pled signal, like c alculating the gravitational p otential for particle 
mesh simulations (Hocknev & Eastwood 1988), or for further pro- 
cessing the discrete galaxy distribution within a FFT based line ar 
Wiener filtering process as described in I Kitaura" &En61in(2008). 

This paper is organized as follows. We begin by briefly outlin- 
ing the general requirements for applying FFTs to real continuous 



signals in section [2] A review of Shannon's sampling theorem and 
the discussion of how to discretize signals sampled to a bounded 
domain are given in section [3] In section [4] we discuss the Fourier 
space discretization of signals, and show, that while it is not possi- 
ble to uniquely recover the entire signal from the sampled dataset, 
it is possible to uniquely restore individual Fourier modes, when 
applying the appropriate filter method. We also show that the dis- 
cretized signal in general will not reflect all physical properties of 
the underlying continuous signal. In section [5TT1 we outline the ideal 
sampling procedure, while in section|6]we discuss the actual prob- 
lems of practical sampling procedures. In section [7] we describe 
our supersampling method, and show how this method can help to 
aliviate aliasing contributions. This is followed by a discussion in 
section[8] 



2 THE REQUIREMENTS OF FFTS 

As already mentioned in the introduction, the FFT is a valuable tool 
in processing huge cosmological datasets due to its computational 
efficiency. This allows us to apply many mathematical operations 
like convolutions or deconvolutions to the data in an computational 
feasible way. However, using FFTs relies on two strong require- 
ments as the function has to be discrete not only in real-space but 
also in Fourier-space. 

It is a well known fact that a natural physical signal like the 
galaxy distribution is neither living in a discretized real nor a dis- 
cretized Fourier-space. The reduction of such a signal to a set of 
finite and discrete sample points, thus, introduces an enormous loss 
of information. 

As will be demonstrated below, the sampling theorem by 
Claude Shannon (Shannori| [l948l.ll949l) requires a function to be 
band limited, meaning, that the function has to be sufficiently 
smooth in order to have its support extending only up to a maxi- 
mal frequency in Fourier-space. If this criterion would be fulfilled, 
then the complete information of the continuous signal can be rep- 
resented by a set of discrete sample points in real-space. However, 
physical theories about the dark matter distribution and the power 
spectrum tell us, that the power spectrum possibly extends over all 
frequencies in Fourier-space, meaning that the matter and galaxy 
distributions cannot be sampled on a discrete grid without loss of 
information. 

In addition, as already mentioned above, discreteness in real- 
space is not enough, and must be understood only as a necessary 
requirement for using FFTs. Using FFTs additionally enforces the 
Fourier-space to be discrete. This is in agreement with the fact that 
the FFT requires the real-space signal to be periodic, meaning the 
signal can be represented by a finite set of Fourier waves. 

For these reasons we claim the requirement of real-space and 
Fourier-space discreteness to be strong criteria in hindsight of in- 
formation conservation, as in general no physical signal will fulfill 
these requirements in a natural way. 

In the following we will discuss how to optimally sample a 
continuous signal in order to conserve as much information as pos- 
sible. In doing so, we will arrive at the ideal instrument response 
function of our computer, which allows us further intuitive insight 
into the problem of representing a real physical quantity on a grid 
of finite sample points. 
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Figure 1. Multiplication a) and convolution b) of a function with the Dirac comb. 



3 DISCRETIZING THE REAL-SPACE 

As already mentioned a necessary requirement for the application 
of FFTs is the real-space discreteness of a signal. However, for the 
case discussed in this work the real physical signal is continuous in 
real-space and hence has to be discretized. This is usually achieved 
by dividing the continuous signal into small regions and associat- 
ing a number with each of those. As this process also involves a 
loss of spatial resolution, and therefore requires to soften sharp fea- 
tures, information gets lost. This process of converting a continuous 
signal into a discrete representation can cause several artifacts like 
aliasing or Gibbs ringing. 

A good way to understand, analyze and eliminate these ar- 
tifacts is through Fourier analysis. In the following we will dis- 
cuss the necessary requirements for real-space discretization, as 
was demonstra t ed by Claude Shannon in his sampling theorem 
(Sh anno nl ll948l.ll949l) . 



3.1 Sampling theorem 

In order to discretize a continuous function, we will represent a 
point sample as a scaled Dirac impulse function. 

With this definition, sampling a signal is equivalent to multi- 
plying it by a grid of impulses, one at each samp le point x,- = /Ax, 
where j is an integer and Ax is the grid spacing dMarschner & Lobbl 
1 19941) . as illustrated in Figure Q] 

Let f(x) be a continuous function, then its sampled version 
g(x) can be expressed by: 



g(x) = Tl(x)f(x), 



(1) 



where the sampling function Tl(x) = 5 D (x — jAx) is a Dirac 

comb or impulse train. With this definition we yield the discrete 
function at the sample positions gj as: 



gj = f(jAx) , 



(2) 



However, in order to demonstrate the origin of aliasing effects 
in the following we will focus on the Fourier transform of equation 



(0. By making use of the convolution theorem [A II which states 
that the product of two functions in real-space yields a convolu- 
tion in Fourier-space, the Fourier transform of the sampled function 
g(p) can be expressed as: 



!(/>) = (n°/)(p), 



(3) 



where the'-symbol denotes the Fourier transform of a function and 
the circle symbol o denotes a convolution. 

As demonstrated in Appendix lA2l the Fourier transform of the 
Dirac comb Tl(p) is again a Dirac comb and is given by: 



JPs) 



(4) 



where p s = 2n/Ax is the repetition length in Fourier-space, which 
is often called the Nyquist rate. Since n(p) is a Dirac comb the 
convolution in equation I0 amounts to duplicating f(p) at every 
point of Tl(p): 

sip) = (n° /)(/>) 



Ps 2 f(P-JPJ 



(5) 



as displayed in FigureQ] This demonstrates that the Fourier repre- 
sentation of the sampled function g{p) is a superposition of the true 
continuous Fourier transform f(p) and all its replicas at positions 
j p s = j(2n)/Ax. We call the copy of f(p) at j = the primary 
spectrum and all other copies alias spectra. This result reveals that 
the sampling operation has left the original input spectrum f(p) 
intact, just re plicating it per iodically in the Fourier domain with a 
spacing of p s JWolberdll997» . 

This suggests to rewrite the sampled spectrum g(p) as a sum 
of two terms, the low-frequency (baseband), and the high frequency 
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(aliasing) components; 



3.2 Low-Pass-Filtering 



g(p) = f(p)+ f<P-jPJ 
= f{p) + F HF (p) 



(6) 



The baseband spectrum is exactly f(p), and the high-frequency 
components, F HF (p), consists of the remaining replicated versions 
of f(p) that co nstitute harmonic versions of the sampled signal 
dWolberdl 19971) . 

A crucial observation in the study of sampled-data systems is 
that exact signal reconstruction from the sampled data requires to 
discard the replicated spectra F HF (p), leaving only f (p), the spec- 
trum of the signal we seek to recover (Wolberg 1997). 

If the different replicas f(p) do not overlap with the baseband 
spectrum, then we can recover f(p) from g(p), by simply multi- 
plying g(p) with a function which is one inside the baseband and 
zero elsewhere, and therefore eliminates the high frequency contri- 
butions, see Figure [2] a). On the other hand, if the high-frequency 
contribution F HF (p) has some overlap with the baseband spectrum, 
there is no way to uniquely recover the original signal f(p) from 
its sampled version g(p), see Figure[2]b). Therefore, the only pro- 
vision for exact sampling is that f{p) m ust be undistorted due to 
the overlap with F HF (p) jWolberdfl997l) . For this to be true, two 
conditions must hold: 

(i) The signal must be bandlimited. This avoids functions f(p) 
with infinite extent that are impossible to replicate without overlap. 

(ii) The sampling frequency p s = (2n)/Ax must be greater than 
twice the maximum frequency p max present in the signal. This can 
be understood by looking at Figure [2] This minimum frequency, 
known as the Nyquist rate, is the minimum distance between the 
spectra copies, each with bandwidth p s . 

The first condition merely ensures that a sufficiently large 
sampling frequency p s exists t hat can be used to separate repli- 
cated spectra from each other l Wol berg||l997h . The second con- 
dition provides an answer to the problem of the sufficiency of data 
samples to exactly recover the continuous input signal. It states that 
the signal can only be recovered exactly when p s > 2 p max , with 
Pmax = PNvaui st being the Nyq uist frequency, not to confuse with the 
Nyquist rate dWolberdl 19971) . 

This is t he sampling theorem a s pioneered by Claude Shannon 
in his papers jShannon|[l948l , l 19491) . 



As seen in the previous section, when a signal is being sampled it 
must be band-limited if we are to recover its information content 
correctly. Natural signals, however, are not generally band-limited, 
and so must be low-pass filtered before they are sampled or equiv- 
alently, the sampling operation must include some form of local 
averaging. 

It is therefore required to apply a filter to the signal, which cuts 
away all Fourier modes higher than a certain maximal frequency 
Pmax = Pj/2. Such a filter is called an ideal low pass filter, and it's 
Fourier representation is given as: 



W(p) = 



for \p\ < p max 
for \p\ > p max 



(7) 



It is ideal in the sense, that it has unity gain in the pass-band re- 
gion, hence not introducing any attenuation of the Fourier modes 
to pass, and that it perfectly suppresses all the power in the stop- 
band region. Applying this ideal low-pass filter to the signal by 
multiplication in Fourier-space yields the low-pass filtered signal: 



h(p) = W(p)f(p), 



(8) 



which according to the convolution theorem IA1I translates to its 
real-space representation as: 



h(x) = (Wof)(x). 



(9) 



The real-space representation W(x) of the ideal low pass filter is 
calculated in Appendix |A3| and yields a sine function: 



W(x) 



■ sincip^x) . 



(10) 



Hence, low-pass filtering is equivalent to convolving with a sine 
function in real-space. 

As the sine function has not such a sharp localization in space 
as a Dirac delta distribution, any sharp feature will be smeared out. 
In this sense low-pass filtering reduces spatial information, which 
cannot be uniquely restored. 

The low-pass filtered function h(x) now meets the require- 
ments of the sampling theorem, and can be discretized in real-space 
with a grid spacing Ax < n/p max . With these requirements the sam- 
pling procedure for a natural signal turns into: 

g(x) = U(x) h(x) = U(x) (W o f) (x) . (11) 

This process is graphically represented in Figure[5]a). 
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Figure 3. Aliased sine function for N=32. The vertical lines limit the base-band. 



3.3 Sampling on a finite real-space domain 

We have demonstrated above, that a band-limited signal can be 
uniquely recovered from its discretized approximation. This is due 
to the fact that there exists a sufficiently high sampling frequency 
to separate the aliased copies from the base-band spectrum and that 
the sampling operation has left the original spectrum f(p) intact. 
However, this is closely related to the fact that the samples are dis- 
tributed over the infinite domain of real-space as the sum in the 
sampling operator IT(x) extends from minus to plus infinity since 
only then the Fourier transform of the sampling operator is again a 
Dirac comb. Unfortunately, this is not the case in real world appli- 
cations, especially since it is generally not possible to evaluate the 
infinite sampling sum with a computer. Therefore, one generally 
wants to restrict the investigation to a finite subset of samples. For 
technical reasons one usually chooses a rectangular sub domain. 
Therefore the sampling operator can be split in a zero centered sub 
domain and the surrounding rest as follows: 

n(x) = Y j 8 D (x-jAx) 



= £ S D (x - jAx) + £ K(* " J' A *) + S °i x + J' A *)] 

= n N (.x) + 5 N (x), (12) 

where we introduced the finite sampling operator TIn(x) and the 
trans-domain sampling operator (Jordan operator) Ejv(x). The sam- 
pled signal on the sub-domain is then given by: 



not be reconstructed uniquely. This can be easily seen by Fourier 
analysis. 

According to Appendix I A4I the Fourier transform of the finite 
sampling operator is given as: 



Tl N (p) = asinc N (p) , 

with the aliased sine function asincN(p) being: 



asincN(p) 



i(jf (N+lj) 



(14) 



(15) 



The aliased sine function is plotted in figure [3] It is obvious that 
the Fourier transform of the finite sampling operator is not a Dirac 
comb anymore, though it still exhibits the property of having res- 
onance peaks at the repetition length p s . Sampling on a finite do- 
main therefore amounts to convolving the Fourier spectrum of the 
true continuous signal with the aliased sine function. This opera- 
tion however does not leave the base-band spectrum intact, and, 
unlike as in the Sampling theorem, unique recovery of the underly- 
ing signal is not possible anymore. It is also obvious that due to the 
shape of the aliased sine function the base-band spectrum will al- 
ways be affected by the higher order spectra, which are overlapping 
with the baseband. However, low pass filtering will not completely 
eliminate these aliasing contributions but it can alleviate them. 

Therefore, in general it is not possible to uniquely restore in- 
formation from signals which have been sampled on finite domains, 
and interpretations drawn from these sampled signals will always 
be afflicted with some uncertainty. 



gjv(x) = g(x) - R N (x) , 



(13) 



with gn{x) = n w (x)/(x) and the rest term Rn(x) = Ejv(x)/(x). As 
demonstrated above g(x) contains all information to uniquely re- 
cover the band-limited signal. Thus, by reducing the set of signal 
samples to a finite subset information gets lost and the signal can- 



4 DISCRETIZING THE FOURIER-SPACE 

In the previous section we discussed the theory of discretizing the 
real-space representation of a continuous signal. We demonstrated 
that in real world applications due to the finite sampling operator 
the sampling theorem is not applicable anymore in a strict sense. 
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However, discretizing the real-space representation of the signal is 
usually not enough for data processing purposes. It can for instance 
be easily demonstrated that the Fourier transform of such a discrete 
signal is still a continuous function in Fourier-space. Nevertheless, 
especially applications in which FFT techniques are used require 
also discrete Fourier-space representations of the signal. Therefore, 
in the following we will derive a sampling method, designed to be 
used with FFTs. 



4.1 FFTs and the Fourier-space representation 

FFTs, since they are based on the Fourier- series, require the func- 
tion to be periodic on the observational interval L. This assumes 
that the observed signal can be decomposed into eigen-modes of 
a resonator with length L. This resonator has a ground frequency 
of po = 2n/L and all other harmonic frequencies can be obtained 
by multiplying with an integer p, = p i. Thus, the Fourier domain 
observation obtained via an FFT is discrete. It is worth mention- 
ing, that this discreteness in Fourier-space arises from applying the 
FFTs and therefore implicitly assuming the function to be periodic 
on the observational interval. In doing so we introduced a second 
discretizing process, but this time the Fourier-space is being dis- 
cretized. Thus, FFTs should be considered as an additional filter 
for the true underlying signal. 

As a consequence the question arises how the Fourier modes 
obtained by an FFT of the observed signal are related to the Fourier 
modes of the true natural signal. This is of special interest in the 
case of cosmological matter power spectrum estimation. To demon- 
strate this we will follow the usual ideal approach to discretize a 
signal, by first low-pass filtering and then sampling it to grid posi- 
tions. The filtered and gridded signal can then be represented by: 



2n/(AxN) = 2n/L and the mode coupling function U{Apk - p): 



8j 



dxW(Ax j-x)f(x), 



(16) 



where W(x) is the ideal low-pass filter as given by equation tlOt . 
Applying the FFT as defined in Appendix |A6| yields: 



2n jk- 



7=0 
N-l 



C ^ e-- n,k -tr I dx W(Ax j - x) f(x) 

A „oo N-l „ce 

a? J dpf " {p) 2 e ~ s "*"* 1 J dx w(Ax j ~ x)e ^~ [ 

|- f~dpf(p) g e- 2 "^ ( e +*P*J) jdS W(x'), 
£.jy p f (p) W(p) g e-^£ (e^»^) 
1 dp f(p) W(P) C £ e 7 "*™ (e V-ipAiyj 



,-^l-ipx 



N-l 



dpf(p)W{ P )cY^(e 

l r 

2n J_ a 



-^fAAxj(Apk-p) 



dpf(p) W(p) U {Apk-p), 



U(Apk-p) = c£ 



,- -fAKxHApk-p) 



N-l 



;=0 

I _ e -V-lAx(A/*-p)A» 



C 



I - e - yf-lAx(Apk-p) 



= Ce 



Ce 



V^T^ (Apk-p) n _ - V^T# (Ap*-p) jv 

sTl!^(Apk-p)(N-\) g e 



e V^T !f (Apk-p) _ e --vCl^S(Aj*-p) 



sin( ^(Apk - p)N) 

(Apk-p) (N-l) U' ' 1 ' > 



sin( A f(A/?yt-p)) 

= Ce-^^^-i'^-^asinciApk-p). (18) 

This immediately demonstrates that the spectral representation ob- 
tained by the FFT is a filtered version of the true Fourier transform 
of the signal. It can be easily seen, that if the true signal is peri- 
odic on the observational domain L then the convolution with the 
filter kernel U(p) vanishes. Hence, in this case we obtain the true 
Fourier representation of the signal. However, as natural signals are 
in general not periodic on the observational domain they cannot be 
decomposed in a discrete set of Fourier waves, and hence some 
Fourier smoothing must be introduce, in order to represent the non- 
periodic function by a periodic function. It is also worth noticing, 
that the mode coupling function does not only affect the amplitude 
of certain modes, but also introduces a phase shifting factor. This is 
due to the fact, that the FFT applies a window which is not centered 
on the origin to limit the real-space domain. For these reasons it is 
not possible to uniquely recover the true Fourier spectrum from the 
FFT Fourier representation of the signal. 



4.2 Sampling in Fourier-space 

In the previous section we demonstrated that applying an FFT or in 
more general a DFT to a low pass filtered non-periodic signal yields 
a Fourier representation, which deviates from the true Fourier rep- 
resentation by a convolution with the mode coupling function U (p). 
The continuous signal therefore has to be filtered in such a way 
that its Fourier representation has a discrete spectrum. This can be 
achieved by sending the signal through a resonator of length L, 
which effectively means the convolution with a Dirac comb in real- 
space. The proof of this can trivially been shown by making use of 
the convolution theorem, and knowing that the Fourier transform of 
the Dirac comb is again a Dirac comb. We introduce the realspace 
replication operator Hr(x) as 



U R (x)= ^ S D (x-jL), 



(19) 



with its Fourier transform, according to Appendix l lA2t . being: 
2/r 



tlR( P ) = — J]6 D (p-jAp). 



(20) 



(17) 



The replication operator allows us therefore to formulate the nec- 
essary step of discretizing the Fourier-space representation of the 
continuous signal, by evaluating the following convolution: 



where we have introduced the Fourier-space lattice interval Ap 



f R (x) = (U R of)(x), 



(21) 
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Figure 4. Sampling process of a galaxy distribution. The galaxy distribution displayed in panel a) is low-pass filtered in panel b) and finally sampled at discrete 
positions in panel c). 



where ,/r(x) is the replicated signal in real-space. The according 
Fourier space representation then reads: 



/*(/>) = n*0)/0) 



In 



^ 6 D (p-jAp)f(p). 



(22) 



This immediately demonstrates that fn(p) is discrete in Fourier- 
space. Forcing a continuous non-periodic function to be periodic, 
or sending it through a resonator, is a filtering process, which dis- 
cretizes the Fourier space representation. We can now replace f(p) 
by the discretized function /r(p) in equation l |17t to yield: 



= dpMp)W(p)U(Apk-p) 



lAp)f(p) W(p) U(Apk-p) 



(23) 



Note, that the mode coupling function now only depends on two 
integers, and therefore is a matrix in discrete Fourier-space. The 
mode coupling function is then: 



Q(Ap(k-l)) = cj] 



N-l 



C 2 



;=0 

= CNdf k , 

where for the last equality we refer to Appendix I A7I 
Using this result in equation {23} then yields: 

h = -rf(kAp)W(kAp). 



(24) 



(25) 



This is a remarkable result. As the ideal low-pass filter W(p), as 
given by equation Q, has unity gain in the pass band, it is possible 
to uniquely recover individual Fourier modes of the true physical 
continuous and non-periodic signal. 

It is therefore clear, that the use of FFT techniques requires to 
sample the function twice, once in real-space, and once in Fourier- 
space. This can simply be achieved by replicating the continuous 
and non-periodic signal and low pass filtering it. According to Ap- 
pendix |A5l the ideal discretization filter kernel ^(x) would then be 
an aliased sine function: 



*P(x) = — asincN(x) . 



(26) 



This filtering procedure can be accomplished in a two step filtering 
process, where one first applies the replication operator and then 
the low-pass filter, or vice versa, to the continuous non-periodic 
signal. From this two step filtering it is also clear that there exist 
two parameters to adjust the sampled field. The first parameter is 
the cut-off frequency p max of the low-pass filter, which controls the 
real-space resolution, the second parameter is the resonator length 
L which controls the Fourier-space resolution, or the spatial domain 
under consideration. 



4.3 The instrument response function of our computer and 
the loss of information 

In the previous sections we demonstrated that in order to process 
information via FFT techniques on a computer the function has to 
be band limited and discrete in Fourier-space meaning the function 
has to be periodic on the observational interval and can be repre- 
sented by only a finite amount of Fourier waves. As natural signals 
in general do not obey these requirements processing the data with 
a computer has to be understood as an additional observational and 
filtering step, which modifies the input data. 

Therefore, in order to make inferences from the information 
stored on a computer towards the information about a real physical 
observable, one must take into account the systematics introduced 
by the sampling process and the use of FFTs or DFTs. The discrete 
representation on a computer is only an approximation to the real 
continuous signal, and hence cannot represent all properties of the 
true original physical observable. This will be demonstrated now 
by the case of the strictly positive galaxy density field. 

In many cosmological applications, like for instance power 
spectrum estimation, galaxies are considered to be point particles 
represented by a Dirac delta distribution. With this definition a 
galaxy distribution consisting of M galaxies can be expressed as: 



p(x) = JV(x-x„), 



(27) 



P =i 



with x p being the position of a galaxy. As can be seen easily this 
density field is strictly positive by definition. However, applying the 
ideal low pass filter yields: 



p'(x) 
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Figure 5. Ideal sampling scheme a) and two stage supersampling scheme b). 



M 

= J^Wix-x,,). (28) 



This can be understood as replacing the Dirac delta distribution by 
the ideal low-pass filter kernel and therefore giving each galaxy the 
shape of a sine function. 

In Figure [4] we display a simple example of a low-pass fil- 
tered galaxy distribution with positions x p e (-8, -3, - 1, 0, 1, 5) to 
be sampled on a grid with grid spacing Ax = 1. As already sug- 
gested above, the low-pass filtering introduces some sort of spatial 
smoothing due to the finite width of the sine kernel. Hence, two 
galaxies being closer together than the grid spacing cannot be re- 
solved independently as can be seen in the cases of the galaxies at 
position x p = - 1 and x p = 0. Another important thing to remark 
is that due to the oscillatory nature of the sine function, the super- 
position of sine functions will lead to positive and negative inter- 
ference. This might introduce density peaks at positions where no 
peak would be observed in the true natural signal, which is demon- 
strated in Figure|4]at position x = -5.5. 

As the sine function is not strictly positive, the low pass fil- 
tered galaxy density given in equation {28} does not possess the 
physical property of being strictly positive as the original observ- 
able, demonstrated as well by Figure[4] This property will only be 
restored in the limit of infinite resolution, when the sine function 
approaches the Dirac delta distribution. 

Being a strictly positive density field thereby is a physical 
property which cannot be represented by the sampled galaxy den- 
sity field. This result will also be true for other physical properties 
which will only be recovered in the limit towards infinite resolu- 
tion. 

However, it is worthwhile mentioning that though the galaxy 
density field may have negative contributions all integral quantities 
like total number of galaxies or total mass, are identical to the ones 
which could be obtained by integration over the true natural signal. 
This is due to the fact that the low-pass filter conserves the zeroth 
Fourier mode of the true natural signal. For this reason it is also 
not possible to fix the negative contributions by cutting them away 
or taking the absolute values of the low-pass filtered signal, as this 
operations will not conserve the zeroth Fourier mode and therefore 
violates conserved quantities like the total number or total mass of 
galaxies. 



5 SAMPLING 3D GALAXY DISTRIBUTIONS 

In the sections above we examined the theory of discretizing con- 
tinues signals and demonstrated some of the problems coming with 
this approach. For simplicity we discussed only the Id case, but all 
results can be extended straight forwardly to the 3d case. In this sec- 
tion we will discuss the practical implementation of sampling meth- 
ods in order to process galaxy distributions via FFT techniques. 

5.1 Ideal sampling procedure 

As already described in section |4~31 a galaxy distribution can be 
expressed as a set of point particles. With this definition a 3d galaxy 
distribution consisting of M galaxies can be expressed as: 

M 

p(x) = Y J S D (x'-x , „), (29) 

where x p is the three dimensional position vector of a galaxy. Note, 
that assuming a galaxy to be a point particle is a strong prior. There- 
fore in more general it is possible to allow for the galaxies to have 
a shape s(x), which describes the physical extend of a galaxy into 
space. The new density field p'(x) is then obtained by simple con- 
volution of equation j29\ with the galaxy shape s(x): 

M 

p'(x) = (s op)(x) = Y J s(5t- xp) , (30) 
p=i 

However, such a modification of the galaxy shape would only be 
necessary if the grid resolution is sufficient enough to resolve in- 
dividual galaxies and their inner structure, which is not the usual 
case for applications as discussed in this work. Therefore, in the 
following a galaxy distribution will always be considered to be a 
set of point particles. 

The galaxy distribution as given by equation l |29l l can be sam- 
pled to a grid according to the ideal sampling scheme displayed in 
Figure|5] which consists of two steps. In the first step the high fre- 
quency contributions of the galaxy distribution are eliminated by 
applying the ideal low-pass filter: 

a; 

g© = (ffop)(f) = ^ff(f-4), (31) 

P =i 

which simply states that the galaxies have been given the artificial 
shape of the 3d sine function W(x). In a following step the low- 
pass filtered distribution will be sampled at discrete positions on 
the grid. For practical implementation these two steps will usually 
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Figure 6. Comparison of the real (a) and Fourier-space (b) representations of the NGP kernel (dashed curve), the CIC kernel (dotted curve) and the TSC kernel 
(dashed dotted curve) to the ideal low-pass filter (solid curve). 



be done in one step, by simply convolving the galaxy distribution 
to the discrete gridpoints: 



Qit = ^ W(A.xFt - x p ) , 



(32) 



P =i 



where Ax is the grid spacing and Ft is the triplet of integers Ft = 
i, j, k. The sum in equation ([32) describes the ideal sampling proce- 
dure, which allows to have the best information conservative rep- 
resentation of the continuous signal in discrete form. However, as 
the ideal low-pass filter kernel extends over all space, the sum over 
all particles has to be evaluated for each individual voxel of the dis- 
crete grid, which makes this procedure impractical for real world 
applications. 



6 PRACTICAL SAMPLING 

Due to the infinite support of the ideal low-pass filter kernel in real- 
space the ideal sampling method (by convolving with such a filter) 
is in general not feasible. For many applications it is computation- 
ally too expensive. For this reason one usually chooses a practi- 
cal approach by approximating the ideal sampling operator, by a 
less accurate, but faster calculable sampling operator. This usually 
means to approximate the low pass filter, by a function with com- 
pact support in real-space. As a result the convolution can be calcu- 
lated faster, for the prize of not completely suppressing the aliasing 
power in the stop-bands. Approximating the ideal sampling opera- 
tor therefore is always a trade off between accuracy and computa- 
tional speed. 

In the literature of Digital Signal Processing this problem 
is well known and studied for many years already. Especially 
the literature of modern 3D computer graphics provides a lot 
of practical solutions to the problem of sampling and filtering. 
Many very detailed studies about the optimal filter approxima- 
tions have been made, and can be found in the signal pro- 



cessing literature dTheuM et alj|200d: iMitchell & Netravali|[l988l ; 
iMarschner & Lobblll 994 IWolberdl 19971) . 

In the following we will discuss some sampling techniques as 
commonly used in cosmology and display their strength and weak- 
nesses. 



6.1 Filter Approximation 

There are sever al known practical issues when dealing with filter 
approximation dMarschner & Lobblll994l) . Using filter or filter ap- 
proximation will in general lead to a variety of effects like smooth- 
ing ("blurring"), which refers to the removal of rapid variations in 
the signal by spatial averaging. Or, one will observe a reshaping 
of the Fourier-space by applying a filter which has not unity gain 
at every mode in Fourier-space. In addition, low-pass filtering step 
discontinuities will r esult in oscillations or ringing, just before or 
after discontinuities dMarschner & Lobdfl994l) . This is the Gibbs 
phenomenon, and is due to the fact, that step discontinuities in the 
signal cannot be represent by a finite superpostion of Fourier waves. 

In cosmology usually t he sub-class of separable lo w-pass filter 
approximations are used ( Hocknev & Eastwoodfl 988). These filter 
obey the relation dMarschner & Lobbl l 1994): 



W(x,y,z) = W s (x)W s (y)W s (z), 



(33) 



where W s (x) is the separated ID filter kernel. This separation al- 
lows for fast com putational speed in performing the 3D convolution 
dFergusonl|200ll) . 

The simplest filter approximation frequently used in cosmol- 
ogy is the Nearest Grid Point (NGP) kernel. This zero- order kernel 
provi des the simplest and fastest interpolation method dDuan et all 
120031) . Each galaxy in the input data is assigned to the value of the 
nearest sample point in the output data. The nearest neighbor kernel 
is defined as: 



1 for -0.5 < x < 0.5 
otherwise 



(34) 



The Fourier response of this kernel is a sine function which 
has a poor localization and pass-band selectivity. This property typ- 
ically leads to low-quality interpolated data with blocking effects 
for signals with high frequency contents like sharp intensity varia- 
tions or high noise level as might be expected in sampling individ- 
ual point particles dDuan et al 120031) . 

Another commonly used filter is the linear filter. This first- 
order kernel linearly interpolates between adjacent points of the 
input data along each dimension dDuan et ai1l2003l) . It is defined 



W S (X) : 



for < \x\ < 1 
for 1 ^ \x\ 



(35) 
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In Cosmology this filter is often referred t o as the Cloud in Cell 
(CIC) scheme dHocknev & Eastw ood 1988). It is popular for pre- 
filtering as it provides a good tradeoff b etween filtering quality 
and computational cost jDuan et al]|2003l) . Nevertheless, a signifi- 
cant amount of spurious stop-band aliasing components continues 
to lea k into the pass-band, contributing to some aliasing JWolberd 
1 1997b . 

The next order kernel is the so called Triangular Shaped Cloud 
(TSC) kernel which is defined as: 



W,(x) 



0.75 - \x\ 2 
0.5(1.5 -W) 2 




for < |x| < 0.5 
for 0.5 < \x\ < 1.5 
for 1.5 < \x\ 



(36) 



dHocknev & Eastwood Il988h . This filter kernel has better stop- 
band behaviour and is therefore usually preferred to the CIC kernel. 

This can be seen in Figure [6] where the Fourier transforms 
of the NGP, CIC and the TSC filter kernel are compared to the 
ideal low-pass filter. In addition to the stop-band leakage of all fil- 
ters they all introduce pass-band attenuation, meaning the power 
of the Fourier modes is suppressed by the filter. This is of special 
relevance in all applications where one is interested in the correct 
spectral representation of the signal, in particular in the case of the 
matter power spectrum estimation. 

This reshaping of the Fourier-space due to the application of 
an imperfect low-pass filter can be undone by deconvolving with 
the Fourier transform of the corresponding filter kernel. This pro- 
cedure in fact returns the correct power at most Fourier modes, as 
can be seen in Figure [7] where we demonstrate the power spec- 
tra obtained after applyi ng these filter approxima tions to a galaxy 
mock catalog created bv lDe Lucia & BlaizoJ ^20071) . The left panel 
displays the power spectrum without deconvolving by the filter ker- 
nel while the right shows the same power spectra but deconvolved 
by the corresponding Fourier transforms of the filter kernels. As 
it can be seen in Figure [7j most of the modes now seem to have 
correct power, and the different filter approximations agree to each 
other at the lower modes. Nevertheless, this procedure will also en- 
hance the aliased contribution leaking into the pass-band region, 
which will severely affect the power spectrum towards the highest 
Fourier modes. In general all Fourier modes of the sampled signal 
will be affected by aliasing. However, natural signal Fourier ampli- 
tudes tend to drop off with higher frequency and therefore the sam- 
pled signal Fourier modes will be affected less and less towards the 
lower frequencies. 

For this reason the power spec trum is usually trusted only up 
to Fourier modes of k c = OJnN/L JCui et al.ll2008n .in cosmology. 
Figure UJ also demonstrates that the NGP method performs poorly 
and should not be considered to be used in accuracy applications. 

Another important thing to remark is that although TSC cou- 
ples 3 3 cells, while CIC only couples only 2 3 , and therefore has a 
broader spatial support, it does not perform much better than the 
CIC scheme. This is due to the fact, that the overall amplitude of 
the ideal low-pass filter drops off in space, and hence the main con- 
tribution of the corrections will come from cells close to the origin. 
As the amplitude of the ideal low-pass filter decreases rapidly it re- 
quires more and more cells to make as strong corrections as were 
obtained by only considering the closer cells to the kernel origin. 
This demonstrates that in order to make further corrections to the 
filter approximation, one requires to couple more and more cells, 
while at the same time corrections will be smaller and smaller. 



7 SUPERS AMPLING 

As already pointed out above, correct sampling is in practice gen- 
erally not feasible, and some approximation must necessarily be 
made. However, for most natural signals, like the galaxy distribu- 
tion, an approximation method can take advantage of the fact that 
signal content generally declines as the frequency increases. This 
observation resulted in the rule of thumb, that t he matter powe r 
spectrum can be trusted at frequencies below k c JCui et al1 l2008). 
The most common anti-aliasing technique in 3D computer graph- 
ics, supersampling, takes advantage of this by sampl ing the sig- 
nal at a frequency higher than the desired sample rate (Gos s & Wvj 
l200d) . A low-pass filter is then applied to the supersamples which 
attenuates or eliminates the frequency content above a threshold so 
that the signal resampled at the target rate exhibits fewer aliasing 
artifacts. 

As supersampling methods are very successful in reducing 
sampling artifacts in 3D computer graphics, we present an adapted 
supersampling method to be used in cosmological applications. 



7.1 Super resolution and downsampling 

The supersampling method consist of two main steps, the super- 
sampling step, in which the signal is sampled at high resolution, 
and the downsampling process, in which the high resolution sam- 
ples are sampled to the target resolution. In our approach we make 
use of FFTs to allow for pass-band attenuation correction and for 
fast and efficient calculation of the overall supersampling method. 
This two stage filtering process is illustrated in Figure [5] and con- 
sists of the following steps: 

(i) supersampling: The continuous signal is sampled to a grid 
with a resolution n times larger than the target resolution. This is 
achieved by applying the CIC or TSC method to allow for fast and 
efficient computation of the high resolution samples. 

(ii) downsampling: The high resolution samples are corrected 
for pass-band attenuation and low-pass filtered. The high resolu- 
tion low-pass filtered samples are then resampled at the target res- 
olution. 

While the supersampling step is in principle identical to the 
method as described in section [6] the downsampling process re- 
quires some more explanation. 

As already described in section [6] using an imperfect filter 
approximation usually leads to pass-band attenuation which has to 
be corrected. 

This is achieved by applying a FFT to the supersampled sig- 
nal, and then dividing by the Fourier transform of the according 
filter approximation. Reducing the supersampled signal to the final 
resolution requires to low-pass filter the supersamples and sample 
them again at lower resolution. 

In our approach low-pass filtering and downsampling is done 
in one step. The low-pass filter step can easily be achieved in 
Fourier-space by discarding all frequencies higher than a given 
threshold. Introducing the ideal low pass-filter in discrete Fourier- 
space as: 



W, 



for-(^-l)<*<*£ 
otherwise 



(37) 
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Figure 7. Power spectra calculated from a galaxy mock catalog created bv foe Lucia & BlaizoJ ('2007) with the three different filter methods NGP (solid line), 
CIC (long-dashed line) and TSC (short-dashed line). Panel a) displays the power spectrum obtained without further processing, while panel b) demonstrates 
the effect of deconvolving with the filter kernel. In panel c) we show the power spectrum taken using the supersampling technique, prefilterd with CIC kernel 
(solid line) and TSC kernel (dashed line) which cannot be distinguished visually. The vertical line indicates the critical Fourier mode k c below which the 
spectrum usually is trusted. 



with Nss being the number of supersampling cells, yields the low- 
pass filtered supersampled signal g ; as: 



8j 



Css J] ^ 
Css ^ fk 



— 



(38) 



/W.V.S . 



where fk is the pass-band attenuation corrected Fourier transform of 
the supersampled signal and Css is the FFT normalization constant 
according to the super resolution Fourier transform. Downsampling 
can now easily be achieved by introducing the number of cells of 
the target resolution N DS = N ss /n and using j = n /, which simply 
results from the fact that the grid spacing of the target resolution 
grid is n times smaller than the super resolution grid spacing. 



Sj 



Css 



E 



r 2nj'k 



(39) 



Sampling the supersampled and low-pass filtered signal to the tar- 
get resolution grid can be done by applying the Rronecker delta 
function 8 k .., to yield the downsampled signal: 



8f 



C s 



= C D 



E 

E 

J"ds 



Infkj 



(40) 



normalization and apply an inverse FFT to obtain the real-space 
representation for the signal. 

The interpretation of this process is easy. As the FFT assumes 
the signal to be periodic and to consist of a finite superposition of 
Fourier waves, the resolution of the supersampled signal can sim- 
ply be reduced by discarding the high frequency waves, resulting in 
a smoother real-space representation of the sampled signal. The re- 
sult of this procedure is a sampled signal at target resolution which 
is greatly cleaned from aliasing effects. 

The ability of this supersampling method is demonstrated with 
Figure [JJ Here we used a super resolution factor n = 2 and calcu- 
lated the power spectrum from the supersampled signal. For the 
supersampling step we applied the CIC and the TSC method. As 
can be seen in Figure [7] the two results for CIC and TSC cannot 
be distinguished. Also it is worth noticing that the aliasing con- 
tribution had been corrected far beyond the usual limit of trust 
k c = Q.lnNIL. 

In Figure[8]we contrast the supersampling method to the TSC 
method by comparing slices through a sampled 3D mock galaxy 
catalogue in a 500 Mpc box which has been sampled regularly at 
128 3 pixels. At first glance, visually there is not much difference 
between the two samples. However, a closer inspection of the high 
density regions reveals that the density field obtained with the pure 
TSC method is more blocky, meaning it can have sudden density 
jumps from one pixel to the next. In the case of the supersam- 
pled density field these high density pixels are surrounded by pixels 
with varying densities. In this way the sharp edges are softened, as 
should be expected by an anti-aliasing technique. Additionally, one 
can observe, that the filaments are better connected in the case of 
the supersampled density field. This is due to the additional super 
resolution information which was taken into account in the super- 
sampling process. 



with C D s being the FFT normalization according to the target res- 
olution. 

Hence, equation ((40} provides an easy prescription for the 
low-pass filtering and downsampling process. Downsampling can 
be achieved by sorting the Fourier modes of the supersampled sig- 
nal up to the Nyquist mode of the target resolution in the target 
resolution grid in Fourier-space, rescale it to the correct Fourier 



8 DISCUSSION AND CONCLUSION 

In this work we reviewed and discussed the use of FFT techniques 
in cosmological signal processing procedures. As described in sec- 
tion[2]the application of FFTs requires the signal to be discrete in 
real and Fourier-space. As natural signals are in general continuous 
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in either spaces, they must be approximated by discretized repre- 
sentations with the constrained of conserving as much information 
of the true natural signal as possible. According to Shannon's the- 
orem this can be achieved by low-pass filtering and sampling the 
continuous signal at discrete positions. 

However, we also demonstrated that Shannon's theorem is a 
necessary but not sufficient requirement, for processing a signal 
via FFTs or in more general DFTs. The application of FFTs or 
DFTs additionally requires the signal to be periodic on a bounded 
spatial domain, or in other words, the signal must be decompos- 
able into a finite amount of Fourier waves. Thus, the application 
of FFT or DFT techniques implicitely assumes the discreteness of 
Fourier-space which is generally not given by natural signals. This 
introduces an additional filtering procedure for the true continuous 
Fourier modes, as shown in section |4~Tl 

In section POl we demonstrated that, due to discretizing, the 
physical property of being a positive density field is not reflected 
by the discretized representation of the continuous galaxy distribu- 
tion. This is due to the convolution with the ideal low-pass filter 
kernel which is a sinus cardinal function and thereby not positive 
definite. The loss of physicality is thus an expression of the loss of 
information due to discretizing the continuous signal, and hence is 
a fundamental problem of the method itself. 

In addition to this fundamental problems, we discussed that 
ideal sampling is in general not applicable to real world problems. 

Signal processing in cosmology or any other field is a techno- 
logical challenging problem governed by the requirement of com- 
putational feasibility of the sampling procedure. As ideal sampling 
is in general computational too expensive to be used for real world 
signal processing applications one usually introduces low-pass fil- 
ter approximations to allow for faster computation for the expense 
of introducing sampling artifacts like aliasing. 



It is, however, often exactly this leakage of aliasing power 
from the stop-band which makes further processing of the sam- 
pled signal difficult. For example estimations of higher-order spec- 
tra with FFT techniqu es will be errone ous due to sampling arti- 
facts. As mentioned in lCui et all J2008t) . so far there is no known 
approach to accurately correct for the sampling effects in measur- 
ing higher-order spectra like the bi-spectrum with FFT techniques. 
Also note, that iterative signal reconstruction processes, which uti- 
lize operations like real-space multiplications, deconvolutions and 
other signal processing operations, may enhance or distribute this 
false aliased power into regions of the power spectrum formerly 
unaffected by aliasing. 

The aim for signal processing technologies therefore is to find 
a low-pass filter approximation which sufficiently suppresses these 
sampling artifacts, while at the same time still being computa- 
tionally far less expensive than the ideal sampling procedure. The 
Digital Signal Pr ocessing litera ture provides plenty of possible ap- 
proaches (see e.g. Smith ( 2002)) which could be introduced to var- 
ious cosmological signal processing problems. 

Recently many authors presented new methods to reduce 
aliasing effects in cosmological po wer spe ctrum estimation due to 
imperfect filter approximation. Jing (2005) proposes to evaluate all 
aliasing sum, and additionally uses an iterative scheme, which re- 
quires a priori knowledge about the true power spectrum , to cor- 
rect for the aliasing effects. More recentlv lCui et al. I d2008l) demon- 
strated that by using scale functions of the Daubechies wavelet 
transformation as an approximation to the ideal low-pass filter 
aliasing can be greately suppressed. This is achieved by allowing 
the filter approximation to extend over a larger support in real- 
space. 

Enlarging the spatial support of the low-pass filter kernel will 
lead to a better approximation of the ideal low-pass filter. Hence, 
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As it tends towards evaluating the sum given in equation < !32b . 
this approach of using filter with greater and greater spatial support 
becomes more and more computationally impractical. 

For these reasons we introduced a new supersampling tech- 
nique frequently applied as an anti-aliasing tech nique in 3D com- 
puter graphics (wolberglll997l:lGoss & WuhoOCh . 

As described in section [3~Tl aliasing is a result of sampling, or 
more specifically, the lack of a sufficient sampling rate. 

Supersampling, as the name suggests, solves the aliasing prob- 
lem by taking more samples than would normally be the case with 
usual particle assignment schemes as CIC or TSC. By taking more 
samples at sub-pixels, we are able to more accurately capture the 
details of the natural continuous signal. The target pixel value is 
then obtained by averaging the values of the sub-pixel reducing the 
aliasing edge effects in the signal. Hence, supersampling reduces 
aliasing by bandlimiting the input signal and exploiting the fact that 
the signal content generally declines as the frequency increases. 

The supersampled signal still contains artifacts due to alias- 
ing, but the artifacts are less prominent as if the signal was sampled 
directly at the target rate. Nevertheless, there are two problems as- 
sociated with supersampling. The first problem is, as already men- 
tioned, that the newly designated Nyquist frequency of th e super- 
resolution samples continues to be fixed I Wolbergl ll997l) . Hence, 
there will always be sufficiently higher frequencies tha t will alias. 
The second problem is cost of memory dWolberg|ll997h . To sample 
a 3D signal at double resolution requires eight times more mem- 
ory. On the otherhand, if memory is no issue, the supersampling 
method is computationally less expensive than introducing a next 
order low-pass filter approximation, while at the same time pro- 
viding better aliasing suppression. In addition, the supersampling 
procedure, as proposed in this paper, incorporates pass-band atten- 
uation corrections, which otherwise would have to be corrected for 
in an additional separate step, when applying just a low-pass filter 
approximation. 

This supersampling technique, therefore, can be understood as 
being complementary to the approach of finding better low-pass fil- 
ter approximations. The combination of better filter approximations 
with the supersampling method will naturally increase the quality 
of the sampled signal further. 

Nevertheless, the problem of signal processing in cosmology 
via FFT techniques is a complex one. Beside being subject to fun- 
damental issues of discritizing physical quantities, it mainly is a 
technological challenging problem, which requires to efficiently 
apply a low-pass filter to an observed continuous signal. The op- 
timal solution for each individual application may also vary from 
case to case, as it can be dependent on the individual properties of 
the underlying signal and scientific application (e.g high accuracy 
power spectrum estimation). 

We therefore believe, that the supersampling method pre- 
sented in this paper will be useful for cosmological signal process- 
ing, permitting highly accurate power spectrum estimates. 



by approximating the ideal sine function closer and closer aliasing 
can be suppressed below any limit. 

This approach is well known in the DSP literature, which 
favors the windowed sine functions as optimal approximations 
to the ideal low pass filter dTheuBl et alJl200Ct iDuan et al.|[2003l ; 



Marschner ~& Lobbll 1 9941; fSmith 1 2002) 
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APPENDIX A: CONTINUOUS FOURIER 
TRANSFORMATION 

We define the Fourier transformation as: 

f(k) = j f(x)e- ik 'dx 

and the according inverse Fourier transform is defined as: 

m = ^J f\ky kA dk 



(Al) 



(A2) 



Al Convolution Theorem 

The convolution theorem states that a convolution in real-space is 
a product in Fourier-space, and a convolution in Fourier-space is a 
product in real-space. This can be shown as follows. 

Let h(x) be the convolution of f(x) with g(x) defined as: 



h(x) = J f(y)g(x-y)dy. 

Applying the Fourier transform lAll to h(x) yields: 

h(k) = J h(x)e- ikx dx 

= j f(y)j e- k >g{x-y)dydx 

= j e- ik > f(y)dyj g(y')dy' 
= f(k)g(k), 



(A3) 



(A4) 



where we substitute with y' = x - y and used lAll to write the last 
line. Hence, a convolution of the functions f(x) and g(x) in real- 
space is a product of the respective Fourier transforms f{k) and 
g(k) in Fourier-space. 

The second part of the prove is analog to the first one presented 
above. 

Let h(k) be the convolution of f(k) with g (k) defined as: 



!<k) = hL f(P'>8(k-p)d P . 

Applying the inverse Fourier transform lA2l to h(k) yields: 



(A5) 



h(x) 



h(k) e ikx dk 



2n X 

l ~ff e ipl f(p)dpj e ik ' x g(k')dk' 



e ikx g(k- p)dpdk 



(2 

fix)g(x), 



(A6) 



where we substituted with k' = k - p and used lA2l to write the last 
line. Hence, a convolution of the functions f(k) and g{k) in Fourier- 
space is a product of the respective inverse Fourier transforms /(x) 
and g(x) in real-space. 



A2 Fourier transform of the sampling operator 

The sampling function n(x) is given as an impulse train: 



m=-oo 

© 2006 RAS, MNRAS 000.mfT7l 



(A7) 



where Ax is the distance between two sampling positions, or the 
sampling interval. Applying the Fourier transform lAll to the sam- 
pling function n(x) yields: 



H(k) 



E 



(A8) 



As can be easily seen, equation lA8l has the form of a Fourier series, 
in fact it is the Fourier series of a Dirac comb except for a constant 
factor, as we are going to show now. The Fo urier series of a period ic 
function f(z) with period Az is defined as iLang & Puc ker 1998): 



m=-oo 

with the series coefficients C m being: 

C„, = -J- f A ~ f(z)e'i""dz. 
Az Jo 



(A9) 



(A10) 



As a Dirac comb n(z) is periodic with period Az, it can be repre- 
sented by a Fourier series as: 



n(z)= 2 S D (z-mAz)= ±- e " ¥l '" Z - 



With this expression we can express fl(A:) as: 
ft(Jfc) = Ps V L e -x^ 



(All) 



(A12) 



where we introduced the periodicity length p s in Fourier-space and 
the uncertainty principle p s Ax = In, Equation |A12] shows that the 
Fourier transform of a Dirac comb is again a Dirac comb. But it 
is important to note that this proof is only valid under the assump- 
tion that the Dirac comb samples at a certain period. As soon as 
the sampling intervals of the Dirac comb become non-periodic this 
proof does not apply any longer. 



A3 Fourier transform pair of the ideal low pass filter 

The ideal low pass filter is defined by: 



W(k) : 



for \k\ < p ma , 
for \k\ > p max 



(A13) 



Applying the inverse Fourier transform l A2I we yield: 



W(x) 



1 f 

27r J-<w 



W(k)e ikx dk 



1 



;27rx 
SiniPmaxX) 



71 X 

Pmax 



sinc{p max x) 



(A14) 



where we have introduced the sinus cardinal sinc(x) = sin(x)jx. 
Hence, the real-space representation of the ideal low pass filter is a 
sine function. 
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A4 Fourier transform of the finite sum sampling operator 

The finite sum sampling operator is given as: 

N/2 

U N (x)= J] S D (x-jAx). (A15) 

j=-N/2 

Applying the Fourier transform lAll then yields: 

N/2 N/2 

tl N (p)= £ e-+*>l*= {e-^"^ ■ (A16) 

j=-N/2 j=-N/2 

Using the closed form of the geometric series: 



a L - a u+ > 



Z« m 



1 -a 



(A17) 



with a = e~ L = -N/2 and V = N/2 yields: 

Ajv(p) 



e ^PT P Ax^ _ g-V^lpA^tf +1) 

J _ g- V-T/>Ax 



sin(^(A^+ 1)) 

sin(^) 
asinc N (p) , 



(A 18) 



where we have defined the aliased sine function asinc N (p) = 
sin ((kAx)/2(N + 1)) / sin ((kAx)/2). 

A5 Ideal discretization kernel 

To discretize the continuous and non-periodic function in real and 
Fourier space, one has to convolve the signal first with the repli- 
cation kernel, and then apply a low-pass filter, or vice versa. Thus, 
the ideal discretization kernel will be a convolution of the ideal 
low-pass filter, with the replication operator. This can be written as 
follows: 



V(x) = j 



d z w(x - z)n R (z) . 



(A19) 



According to the convolution theorem we can express the Fourier 
transform ^(p) as: 



f(p) = w{p)n R { P ) 

2n 

T 



Applying an inverse Fourier transform then yields: 

i z X e ' p * sD ip ~ j Ap) ^ (p)dp 

j=-co 
, M/2 

t z (^r . 



(A20) 



(A21) 



j=-M/2 



where we made use of the fact, that the ideal low-pass filter given 
by |A13| vanishes outside the base-band. One can now use the closed 
form of the geometric series given in equation JA19t to yield: 



j -ixApM/2 _ ixhp(M/2+\) 

L 1 - e ixh P 



1 e' x ^ e -i^P(M/2+\) _ g«A*(M/2+^) 



£ t — e 



I sin(x^(M+ 1)) 
—asinc M (x) . 



(A22) 



Thus the ideal discretization kernel is an aliased sine function in 
real-space. 

A6 Discrete Fourier transformation 

we define the numerical Fourier transformation as: 

N-l 

j^C^x/ 2 "^ (A23) 
j=o 

and the inverse Fourier transform 



i N 
x ] = cY aW ?"*& =C J w-' 



(A24) 



The normalization coefficients C and C are chosen such that they 
fulfill the requirement: 



cc=I. 

N 



A7 Discrete mode coupling function 

Here we proof the equality 



(A25) 



CNSf,, 



(A26) 



by evaluating the sum. First we discuss the case j ■$■ i where we 
will consider the sum 



Z< 



2nk^-(j-i) 



Z2n I — Ln 
cos —k(j - i) + V- 1 sin —k(j - i) 



In. 



N 

N-l 



Z2ti i — v-i 2tt 
COS —k(j - i) + V-l 2j sin -77*0' - 



A- 1 



(A27) 



By making use of the trigonometric sums: 

^ sin((iV-i)x) j 

^cos(fa)= iL-^ >- -- 



and 



2sin(f) 



A "' sin(^x)sin(fx) 



£j sin(foc) : 



we yield: 



sin(f) 



(A28) 



(A29) 



l + Yj cos -77*0 - + "^-i Z sin 17 k( -! ~~ !) 



sin 



1 + 



1 



2sin(§(y-i)) 2 
sin - 0) sinWj - (')) 



(A30) 
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because of sin(n(j - j)) = for all j and i the last term vanishes, 
and we yield 

1 sin(2^0-0-f0-0) 

= o + \ ■ ^ • (A31) 

2 2sin(f0-i)) 

Since the sinus is 2^-periodic we can write: 

sin - - £ (J " ')) = - sin ( £ ( j - /)) (A32) 

Using equation [A32]in equation lA31l we yield: 

1 sin(fQ-Q) 

2 2sin(f(y-i)) 

= |-|=0- (A33) 

Hence, for the case j ^ i the sum will always be zero. Now we 
discuss the case j = i, where the sum yields: 

N-l N-l 

je»*^W = i]l = JV, (A34) 
Hence, we end up with the result: 

s £!;{ <**> 

Which yields: 

J «"* H = MSg . (A36) 
Thus, equation l lA26t has been proven. 

This paper has been typeset from a Tj5X/ KTgX file prepared by the 
author. 
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